{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0c1ed389-03b2-4725-a908-03180a76c578",
   "metadata": {},
   "outputs": [],
   "source": [
    "suppressMessages(library(ArchR))\n",
    "suppressMessages(library(Seurat))\n",
    "suppressMessages(library(Signac))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fe2d3533-2c68-42df-99d0-aaaf77544d61",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Successfully loaded ArchRProject!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "proj <- loadArchRProject(\"../data/VisiumHeart\", showLogo = FALSE)\n",
    "metadata <- as.data.frame(proj@cellColData)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5e354925-a7ee-42ec-ab7d-98593a6dd847",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>319939</li><li>17700</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 319939\n",
       "\\item 17700\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 319939\n",
       "2. 17700\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 319939  17700"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "counts <- readRDS(\"../data/VisiumHeart/PeakMatrix.Rds\")\n",
    "dim(counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6ba803bd-41eb-4d42-a039-762fb553be58",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "An object of class Seurat \n",
       "319834 features across 17700 samples within 1 assay \n",
       "Active assay: peaks (319834 features, 0 variable features)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "chrom_assay <- CreateChromatinAssay(\n",
    "  counts = counts,\n",
    "  sep = c(\"_\", \"_\"),\n",
    "  genome = 'hg38',\n",
    "  min.cells = 30\n",
    ")\n",
    "\n",
    "obj.atac <- CreateSeuratObject(\n",
    "  counts = chrom_assay,\n",
    "  assay = \"peaks\",\n",
    "  meta.data = metadata,\n",
    "  names.field = 1, \n",
    "  names.delim = \"#\")\n",
    "\n",
    "obj.atac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8f52ded9-f372-4161-b7d6-b798c4c50851",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Harmony 1/10\n",
      "\n",
      "Harmony 2/10\n",
      "\n",
      "Harmony 3/10\n",
      "\n",
      "Harmony converged after 3 iterations\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(reticulate)\n",
    "sc <- import(\"scopen\")\n",
    "matDR <- sc$Main$scopen_dr(counts = obj.atac@assays$peaks@counts,\n",
    "                           verbose = 1)\n",
    "\n",
    "matDR <- t(as.matrix(matDR))\n",
    "\n",
    "colnames(matDR) <- paste0(\"PC_\", 1:ncol(matDR))\n",
    "rownames(matDR) <- colnames(obj.atac@assays$peaks@counts)\n",
    "obj.atac@reductions[['scopen']] <- CreateDimReducObject(embeddings = matDR,\n",
    "                                                        assay = \"peaks\",\n",
    "                                                       key = \"PC_\")\n",
    "DepthCor(obj.atac, reduction = \"scopen\", n = 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "293182e8-da11-41df-9445-797745d11421",
   "metadata": {},
   "outputs": [],
   "source": [
    "obj.atac <- RunUMAP(object = obj.atac, \n",
    "                    reduction = 'scopen', \n",
    "                    dims = 1:30, \n",
    "                    verbose = FALSE)\n",
    "\n",
    "suppressMessages(library(harmony))\n",
    "\n",
    "obj.atac <- RunHarmony(\n",
    "  object = obj.atac,\n",
    "  group.by.vars = \"Sample\",\n",
    "  reduction = 'scopen',\n",
    "  project.dim = FALSE,\n",
    "  assay.use = 'peaks',\n",
    "  plot_convergence = FALSE,\n",
    "  verbose = TRUE\n",
    "    \n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a3e212a-36be-4d92-972c-84447ebe892b",
   "metadata": {},
   "outputs": [],
   "source": [
    "obj.atac <- RunUMAP(obj.atac, \n",
    "               dims = 1:30, \n",
    "               reduction = 'harmony',\n",
    "               reduction.name = \"umap_harmony\",\n",
    "               reduction.ke = 'umapharmony_',\n",
    "              verbose = FALSE,\n",
    "                   min.dist = 0.4)\n",
    "\n",
    "options(repr.plot.width = 12, repr.plot.height = 5)\n",
    "\n",
    "p1 <- DimPlot(object = obj.atac, reduction = \"umap\",\n",
    "              group.by = \"Sample\") +\n",
    "    xlab(\"UMAP1\") + ylab(\"UMAP2\")\n",
    "\n",
    "p2 <- DimPlot(object = obj.atac, reduction = \"umap_harmony\",\n",
    "              group.by = \"Sample\", shuffle = TRUE) +\n",
    "    xlab(\"UMAP1\") + ylab(\"UMAP2\")\n",
    "\n",
    "p1 + p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2522888b-0193-4f44-babb-13f25eed13b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "saveRDS(obj.atac, file = \"../data/VisiumHeart/snATAC.Rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c63f93cf-e9aa-4c59-b24a-4cbedf1fb411",
   "metadata": {},
   "outputs": [],
   "source": [
    "sessionInfo()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
